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1 . INTRODUCTION 


The key property of elliptic systems is that their solution tends to be 
as smooth as the data and other factors permit. This is in striking contrast 
to hyperbolic systems where singular behavior (e.g., shocks) can arise even if 
all inputs are smooth. 

To illustrate this, consider the following model problem defined in a 
bounded region 0 of if 1 : 

(1.1) div[a grad<}>] - f in 0 

(1.2) b[<j> ] = g on an. 

The boundary operator B has the form 

(1-3) B [<J> ] = a 

where v_ is the outer normal to 9ft. Thus the inputs to this system are 
the coefficients a, a, 3, the data f, g, and the region ft. If all of 
these are smooth, then the same is true of the solution <J> • 

This means that singular behavior can arise only in the cases where the 
inputs are irregular. The following are examples of technical importance. 

Case 1 : Irregular Boundaries . The study of fracture and crack 

propagation involves elliptic systems of various orders depending on the type 
of problem being modelled [1] - [2]. Torsion problems are second order (like 
(1.1)), while plates involve fourth order equations, and shell problems are 
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even of higher order [3] - [4]. Nevertheless, the region ft in question 
generally contains a slit as in Figure 1.1 with the point P being the crack 
tip* In the linear, second order case (a = 1 in (1.1)) it can be shown that 
with Uirichlet boundary condition (a = 0, 0 = 1 in (1.2)) the solution <j> 
behaves like 

(1.4) $ = or^^ sin(0/2) + o(r) 

near the crack tip (r+0). Thus the gradient is singular at the crack tip: 

(1.5) |v<$> | = 0(r ) as r*0. 

The coefficient a of the singularity is of great technical importance, and 
is called the stress intensity factor. Knowing it permits estimates on crack 
behavior via energy release rates [5] - [6]. 

While singular behavior at corners in solid mechanics is of great 
technical significance, in fluid mechanics corner singularities are often 
irrelevant artifacts. For example, linear potential flow over a flat plate 
involves the same equations and geometry as discussed above [7]. In this 
case, the gradient V<J> is the velocity field and hence (1.5) predicts a 
square root singularity at P. This does not occur in real flows and is an 
artifact of the linearization. A proper nonlinear analysis shows the velocity 
singularity is a good deal milder. This and other nonlinear effects are 
discussed in Section 3. 
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Figure 1.1. The Slit Region ft: (r,6) are polar coordinates at P 


Case 2: Discontinuous Boundary Operators. Singularities can arise when 
the boundary conditions (1.2) abruptly change type, or what is the same, when 
the coefficients a, g in (1*3) are discontinuous. An example is given in 
Figure 1.2. Observe that in this case the Neuman condition 



is a symmetry condition, and hence this case is exactly the one represented in 
Figure 1.1. In particular, the singularity is described by (1.4). Most of 
the applications in fluid and solid mechanics involving discontinuous boundary 
operators arise in this manner. 



Figure 1.2. Discontinuous Boundary Conditions 
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Case 3: Discontinuous Coefficients * Discontinuities in the coef- 

ficient a in (1.1) can also generate singularities in the solution <J> . 
Cases of technical importance include diffusion problems in regions ft con- 
sisting of different materials [8] - [9]. A typical example is shown in 
Figure 1.3. Using polar coordinates (r,0) at the point P, it can be shown 
that the singularity behaves like 

(1.6) <(>=or^$(0) 


plus higher order terms in r for a suitable function $. The 
depends on the (constant) values aj of a in each of the regions 
ftj (j - !,•••, 4). In general, these singularities can be far more 
than the square root singularity in (1.4). In particular, X > 0 
made arbitrarily small for appropriate choices of a^,.«.,a^ [9]. 


exponent 

serious 
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a 3 

P 
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Figure 1.3. The region ft for the apriori bound (2.22) 


Case 4’: Nonsmooth Data . Cases where for g in (1.1) - (1.2) are not 
smooth arise in solid mechanics when bodies are subject to random loads 
[10]. To date, these have been best treated using the techniques of stochas- 
tic differential equations since singular behavior of the solution <j> tends 
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to be distributed over the entire region ft. These techniques are not 
treated in this paper, and the reader is referred to [11] for more details. 

The presence of singularities, such as described above, tend to 
significantly affect the rates of convergence for both finite difference as 
well as finite element schemes [12]. Except for Case 4 cited above, two 
approaches have been used to treat the singular behavior. One, generally 
called the singular element method because of the way it is used in finite 
element schemes, attempts to incorporate the singular behavior into the 
approximation. The other approach uses grid refinement. 

To put this paper into proper perspective, it should be clearly noted 
that, in general, grid refinement is the best approach. This is particularly 
the case when adaptive strategies can be incorporated into the approximation 
[13] - [15]. Singular elements tend to be useful only in very special 
circumstances, some of which are cited below. However, the methodology used 
to derive singular elements still remains quite relevant. The point is a 
priori knowledge of the nature of the singularity can be of practical benefit 
even if this information is used only indirectly. This is certainly true of 
grid refinement as well as h, p, and h-p versions of the finite element 
method [15] - [16]. Because of this, these methodologies will receive the 

primary emphasis in this paper. 

Sections 2 and 3 are devoted to the types of singular behavior that can 
arise in elliptic systems. Section 3 considers nonlinear effects, and this 
material is apparently new. The final section concentrates on practical 
issues associated with singular elements along with selected numerical 


results. 



2. SINGULAR BEHAVIOR IN LINEAR SYSTEMS 


For simplicity attention will be confined to problems with corner 
singularities. Interface problems (Case 3 in Section 1) are treated with 
similar techniques. 

In the linear case the starting point is the construction of explicit 
solutions using separation of variables for constant coefficient problems and 
simple geometries. Then using well developed techniques from the theory of 
partial differential equation (e.g., modifiers and frozen coefficients 
iterations) one can analyze problems with variable coefficients and rather 
broad classes of regions ft ([17] - [18]). 

Apparently, the first researcher to realize that the form of 
singularities can be obtained by a local separation of variables was L. 
Williams [19]. To describe the results in this classic paper, consider the 
sectorial region shown in Figure 2.1 letting (r,6) denote polar coordinates, 
then (1.1) (with a = 1) becomes 


( 2 . 1 ) 


a 2 * l a* i a 2 * _ 

7 ? 7S7 * 


Assuming a solution of the form 


( 2 . 2 ) 


♦ - r A $(0), 


we are led to an eigenvalue problem for the exponent X and the function 
*: 



X 2 * = 


0 . 


(2.3) 
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The type of singularity obtained thus depends on the boundary conditions. 
Dirichlet conditions 


(2.4) *(0) = *(0 O ) = 0 
yield 

(2.5) r* a sino0, sina0Q = 0. 


The relevant solution is the one with the smallest positive index a. This 
gives 


( 2 . 6 ) 


IT 

. ,1T0* 

= r sin(s-) 

e o 


as the dominant singular term. Observe that V<J> is finite at the corner 
point P, if and only if 0 < 9^ < v; i.e., reentrant corners (0^ > tt) 

yield unbounded gradients. In the case of a crack shown in Figure 1.1, we 
have 0 q = 2ir , and hence (2.6) reduces to the square root singularity given 
by (1.4). 

It is important to note that the order of the differential operator 
exerts an important influence on the type of singularities that are 
obtained. For example, consider the fourth order equation 


(2.7) 



defined in the section shown in Figure 2.1. Assuming a solution of the form 
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( 2 . 8 ) 


Ip = r X+ 1 *(e), 


we obtain 


(2.9) *(0) = a^sin(X + 1)0 + a 2 COs(X + 1)0 + b^sin(X - 1)0 + b 2 Cos(X - 1)0. 


The constants aj , bj depend on the boundary conditions to be imposed on 
0=0 and 0 = 0^. For example, along a clamped radial edge one has 


( 2 . 10 ) 


- n 1 du> n 

“ °- Th °’ 


while a simply supported edge gives 


( 2 . 11 ) 


m = 0, c 0 = 0, 


where 


( 2 . 12 ) 


i a 2 
1 3 to 

* 77 ~ 


2 

1 3u) 3 a) 


v being the Poisson ratio. Substituting into either of the boundary 
conditions gives a nonlinear equation for the exponent X. Except for 
special cases this equation must be solved numerically, i.e., explicit 
formulas for the exponents are not known. Nevertheless, a number of 
qualitative features of solutions to these equations are known (along with 
specific numerical values in technically important cases [19]). 
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An important point is that solutions A can be complex. This implies 
that oscillatory behavior in the radial directions can occur. This is 
particularly relevant for cases where the stress intensity factors are 
approximated using only nodal values of the solution <}> . 

In the case of a crack (Figure 1.2) it can be shown that A = . Thus 
the stress o Q given by (2.12), which is the analog of the gradient for the 
fourth order case, displays a square root singularity. 



Figure 2.1. The sector ft 

Another issue of technical importance concerns the effect of dimension. 
The direct generalization of William's work to three dimensions involves the 
conical region shown in Figure 2.2. Letting (r,0,if /) denote spherical 

coordinates, -one seeks a solution of the form 

<|> = r A $(e,i|0. 


(2.13) 
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In the second order case, one is again lead to an eigenvalue problem except 
now for the Laplace-Beltrami operator on the boundary surface of £1: 


(2.14) 


2 

1 9 4 > 

2 2 
sin 6 9 


Hsr If (sl " e w> + 1(1 + 1)4 ‘ “• 


A further separation of variables gives 


(2.15) 


$ (6 ,ip ) = sin(mij) + ct)P(0), 


where P satisfies the Legendre equation 


(2.16) 



- y 2 ) ^0 + EX (A + 1) - 



0 


with y = cos8 . In the case of Dirichlet boundary conditions, solutions are 
obtained by requiring 


(2.17) p (u 0 ) “ °» Uq = cose 0 . 

Properties of these solutions have been studied for the case 0 < 6 q < tt 
[ 20]. The dominant singularity is independent of ij>(i.e.,m = 0, a=y in 
(2.15)). In this case (2.16) can be solved with Legendre functions P * P^ 
and (2.17) reduces to a nonlinear equation for X. Many of the character- 
istics of the two dimensional case reappear here. For example, |V<f>| is 

finite at P only in the case of convex region (0 < 9 q < y) . 
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Interestingly , this analysis runs into trouble in the case of a crack 
(*) 

0q = tt. Here, (2.17) does not have solutions with m = 0 since the only 

solutions to (2.16) with m = 0 which are finite on -1 u 1 are 

Legendre polynomials P^U “ n ) which satisfy | (±1 ) | = 1 . The net 

result is the question of the completeness of the functions in which the 
singular behavior is described. Nevertheless, if one were to accept the 
limiting behavior as 0 +tt as valid, then the apparently ubiquitous square 
root singularity reemerges. 




Figure 2.2. Spherical coordinates and the conical region ft. 


The case 8q = tt is exactly the one where the region w fails to satisfy an 
exterior cone condition [17] - [18]. Most of the results from the theory of 
partial differential equations as well as embedding theorems for Sobolev 
spaces are not valid in the absense of this condition. In two dimensions, 
this appears to be only a technical point limiting the analytical techniques 
used, while in three dimensions it seems to be a fundamental issue. 



! 


- 12 - 


The open question of completeness of the functions used to describe 
singular behavior in three dimensions is also present in regions other than 
cones • Thus a rigorous account of how solutions behave near planar slices 
in ft, a case of considerable technical importance is still not completely 
known • 

Elementary solutions such as those described above can be used to obtain 
results for rather general regions ft and for problems with variable 
coefficients* Typically this approach will not yield explicit behavior, but 
rather shows that the solutions lie in appropriate Sobolev spaces. With most 
grid refinement techniques, this information is almost as useful as knowing 
the explicit behavior. 

To fix ideas, consider the second order case (1.1) - (1*2), where a is 
a smooth function of the spatial variables and a = 0, 6 * 1 (Dirichlet 

boundary conditions). Let ft be the region shown in Figur 2.3 which has a 
corner at P but is otherwise smooth. Note that the edges leading to P do 
not necessarily have to be straight. For this region we introduce the 
following weighted Sobolev spaces in terms of the radial distance r to P. 
Indeed, let 

(2.18) Ikon = {/ r 6 |oj 1 2 } 1/2 

U,P S 

be a weighted L 2 for the given exponent 8 £ R The higher order spaces 
involve derivatives. Thus 

B - {/ r e (M 2 +( 0 2 ]} 1/2 

1,B a 


(2.19) 


II 0) 
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with analogous definitions for llwl! (t > 1). A priori bounds on the 

t >p 

solution were first obtained in a fundamental paper by Kondrat'ev [17] 

(see also [18]). A typical result in the case of homogeneous boundary 
conditions (g = 0 in (1.2)) shows that if 

t + i - (it /e 0 ) < B < t + 1 

and 

( 2 . 21 ) 11 f 11 1 ,8 < ’ 

j 

then there is a constant 0 < C < », depending only on t,8,fi, and a, 

! 

for which 

( 2 . 22 ) 

The analysis requires that ft satisfy an exterior cone condition and hence 
the interior angle 0^ is restricted to 0 < 0^ < 2 tt . This inequality 

I can be verified directly in the cases where the explicit singular behavior is 

known . 

I 

t 



l * l t+ 2 ,e - clfl t,B 


( 2 . 20 ) 


Figure 2.3 
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The practical use of a priori bounds like (2.22), and well as use of the 
explicit singular behavior (when known) will be discussed in the last section. 


3. SINGULAR BEHAVIOR IN NONLINEAR SYSTEMS 

In this section, we shall study the effects of nonlinearities on singular 
behavior. In this setting, the basic paradigm underlying the analysis in the 
linear case, namely separation of variables, is not available, hence different 
techniques are needed. 

To fix ideas, we consider (1.1) with Dirichlet boundary conditions 
(a « 0, B = 1) in the sector shown in Figure 2.1. We assume the coefficient 
depends on the gradient V<j> as 

(3.1) a = a( | V<j) | 2 ) . 

For the problem to be elliptic, we need 

(3.2) ^ [a(t) 2 t] > 0 
uniformly throughout 0. 

An early approach to this problem was given by Tolksdorf [21]. He looked 
for solutions of the form 

(3.3) + = r X $(6), 

and showed that X , $ are related through a nonlinear eigenvalue problem. 
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A key parameter in the latter is the number 

(3.4) q = lim [t ^ll/ a (t)]. 

As will be seen below, this number determines the nonlinear effects on the 
singularity. Observe that q = 0 in the linear case, and for this system to 
be elliptic it is necessary that 

(3.5) q > ~ • 


Because of the nonlinearities, it is not possible to get explicit 
solutions to the nonlinear eigenvalue problem. However, using maximum 
principle arguments with comparison functions, Tolksdorf was able to obtain 
the following bounds on the minimal exponent X in (3.3): 

(3.6) VV - * - X 2 (0) ’ 


where 


(3.7) 


Xj(e) 


qa Q + / (2q+l)(4q+l) 
a Q (4q + 1) 


and 


(3.8) 


X 2 (0) 


qa, 




(2q + 1) + 


2 2 


a 0 (2q 


+ 1 ) 


with 6 q = 
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Unf ortunately , as will be shown below, the inequalities in (3.6) are not 
strict except in the linear case q - 0. Nevertheless, the results do contain 
useful information. For example, 0 < X < 1 for reentrant corners 
(tt < < 2tt), and hence | V<f> | is singular in these cases. Moreover, this 

is the only property that survives from the linear case, and the nonlinear 
terms definitely affect the asymptotic behavior as r*0. It is also 
noteworthy that this approach can also be used in three and higher dimensions. 

Sharp results can be obtained for the case of two dimensions. The 
starting point is a fundamental paper by L. Lehman [22]. This work was 
limited to the linear case, however, the function theoretic approach used did 
not require the explicit construction of solutions via separation of 
variables. The exact form of the singularity was obtained, but in an indirect 
manner using appropriate transformations and analytic continuation. This is 
why the results are limited to planar regions; however, they do point the way 
to an approach for the nonlinear case via a hodograph transformation. 

In particular, we let u,v denote the components of the gradient V<j> 
and rewrite (1.1) as 


(3.9) 


(a + 


aV> It + <2a,uv) It + (2aV “ v) It + (a + aV) n 


(3.10) 


3u 3v 
37 " 37 


0 , 


where 


V da 
3 = dt 


The hodograph transformation takes the form 


(3.11) 


x = x(u,v). 


y = y(u,v) 
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At points where the Jacobian of the transformation 


(3.12) 


J 


3 u 3 v 3u 3 v _ 3 ^<f > 9 ^<J> 9 

^3y~7y33r = — ITT' ttSy 

9x 9y 


is nonzero, we have 

(3.13) (a + 2a? u 2 ) - (2a^uv - (2a^uv) + (a + 2a? v 2 ) ~ = 0 

«•“> if-SJ-o- 


Because the last relation, a potential 


(3.15) 



i|) can be introduced so that 



This gives 


(3.16) 


V 2 

(a + 2a u) 


d 2 ip 

9 v 2 


, v d 2 ip 

- 4a uv 


+ , + , v 2. a 2 t 
Sh t(st2,,> r2 

d U 


Thus the hodograph takes the original nonlinear problem into a linear one. 
Observe that ellipticity of the equation is not affected by the 
transformation. 

Hodograph transformations have been used in a number of problems [7]. 
The technical difficulty has typically been in deriving boundary conditions in 
the hodograph plane. In fact, the traditional approach in fluid dynamics is 
to let the boundary conditions so obtained actually define the problem that is 
being solved [7, Ch. XX]. Here no such problems arise, as can be anticipated 
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from a careful study of Lehman's work. The cases of a reentrant (0^ > tt ) or 
convex (0^ < tt ) corner are different, so for brevity we consider only the 
former since it is the most important. 

The boundary transformations under the hodograph are given in Figure 
3.1. In particular, the sector ft: 0 <_ 0 < 0^ gets mapped into a double 
cone ft with angle 6^ - tt . 



Figure 3.1. The hodograph 


Introducing polar coordinates 


(3.17) u = pcosu), v * psinoj. 


It can be shown that (3.16) reduces to 


(3.18) ( 2_ ) + I - 0. 

a + 2a V p2 ° 7 7 ? 

* 

Moreover, Dirichlet boundary conditions \j> = 0 on 3ft can be used. 
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Tolksdorf's work shows that V<J> is singular at the corner point r = 

0. Thus in the hodograph plane it is the behavior at p = ® which is 
relevant. Since the coefficient in (3.18) is smooth, Kondret'ev's work [17] 
shows one can freeze it in the neighbor of the singularity; i.e., replace it 
with 


(3.19) 


lim (- 
t+°o a + 2a’ p 


' , V 2 } “ 1 + 2q » 
+ 2a d 


where q is the index (3.4) appearing in Tolksdorf's work. The resulting 
problem 


(3.20) 


‘mi 1 J? 17 


can now be solved by separation of variables. In particular, letting 


(3.21) 


. a . , unr v 

= P sxn(^ -) 

e 0 " ir 


we obtain 

(3.22) a = ± a n , a n = q + / q 2 + — - n , 

00 (a 0 - i r 

where 0^ = a^ir. Since it is the behavior at ® which is relevant, we 
take a = Expressing the results in terms of the radial coordinate 

r in the physical plane, we have 


(3.23) 


|Vt|> | ~ r 
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Hence 


(3.24) 


V<fi~r 


° 0 /(a 0 + l) ~ 1 




or what is the same the exact exponent X in (3.3) is given by 


(3.25) 


X = 


1 + °0 ’ 



+ 



^ 2q + 1 

+ 2 • 

<“o - ” 


Observe that in terms of the lower bound X^ and upper bound X 2 cited 
above we have 


(3.26) Xj < X < X 2 

except in the linear case (q = 0) where 


X 


1 


= X 



Consider now the case of a slit where « 2. The exponent \ reduces 

to 


(3.26) 

Observe that if q is negative, then 


X < 


1 

1 


I 


(3.27) 
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i.e., the singularity is stronger than that occurring in the associated linear 
problem. The limiting case is q = - 2- , and the problem is no longer 

elliptic for smaller values of q. That this type of analysis fails is to be 
anticipated. Indeed, the equations will be hyperbolic in such cases, and the 
presence of real charcteristics means that the type of singularity obtained 
depends on information in the far field away from the singular point. 

If, on the other hand, q > 0 we have a situation where the nonlinear 
terms smooth the singularity; i.e., 

(3.28) 1 ) y • 

A case of interest is compressible potential flow [7] where 

(3.29) a(t) 6 ” 1 = M^[Cq - ( g -~ — L)t]. 

It is readily seen that 

(3.30) q - 
in this case, and so 

(3.31) X = i-Jji . 

For dry air 6 = 1.405 and X is approximately .9. As far as numerical 
approximations are concerned, this is a far less serious singularity than 
X * .5 predicted for linear flows. 
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4. PRACTICAL CONSIDERATIONS 

To fix ideas, consider the planar region shown in Figure 2.3, and the 
linear second order case. 

As far as (nonadaptive ) mesh refinement is concerned, the starting point 
is the identification of a region about the singular point where mesh 

refinement is used. An important practical consideration is that with local 
mesh refinement the number of points used in need only grow like 

n 

0(|lnh|j, where h is the far field mesh spacing. The key to this is the 
existence of an a priori bound like (2.22). For simplicity, we take the 

case t ■* 0 so 

(4.1) 1 - ir/6 0 < 6 < 1. 

The grid goes as follows. For triangles adjacent to P the diameter 6, 
satisfies 

1 

(4.2) 5 ( constant)!^ ~~ ® . 

This is increased by a constant factor until the far field uniform grid with 
spacing h is obtained. Grids of this type have been used in [23] for 
standard finite element methods and in [4] for those based on least square 
approximations . 

It can be shown that for such grids the approximation $ using linear 
elements satisfies 


(4.3) 


{/!*(♦ - «(> n )i 2 } 1/2 < C h {/r 2 6 |D 2 <{> 1} 1/2 
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Numerical evidence shows that if approximate difference quotients are used, 
then approximations to the stress intensity factor a can be 

obtained satisfying 

(4.4) a - a. = 0(h 3/2 ) 

n 

in the case of a slit region (Figure 1.1) [25]. It is a point of practical 

significance that the form of the singularity must be known in order to 

establish the relevant difference approximation for a . Path independent 

n 

J-integrals [26] can be used as an alternative, however, they produce 
approximations which converge at lower rates, typically 0(h) when using 
linear elements. 

In a singular element approach, one seeks approximations in the form 


(4.5) 





where <J> ^ is the (known) singular function, is the intensity of the 

singularity (to be computed), and Nj are th e piecewise linear nodal 

functions* Since variational methods typically yield best approximations (in 

suitable norms), the term a,<b in effect substracts out the singularity; 

h s 

i.e*, for a suitable number or (the exact intensity) 


(4*6) 


*♦ " a V2,0 < *’ 


and thus . 4> “ o4» g can be approximated by piecewise linear functions on a 
regular grid. Stated differently 
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(4.7) 

This means the 


inf .* - n+ s - I E 1 N i l 1 0 < ch I* -M,l 2 0 . 

4^ computed by a finite element approximation will satisfy 


(4.8) 


- V.,o < Ch - 


In addition, the intensity in (4.5) has been observed to converge at 

the rate found in (4.4) [25]. 

Substitution of (4.5) into a variational principle typically leads to a 
matrix problem of the form 


(4.9) 


1 

1— 1 

K , ~ 
Is 


r 

-©-1 

i 

it 

~f~ 

K , 
si 

1 

W 

CO 


b 

i 


f 

s_ 


where is a regular finite element matrix, K^ g contains inner products 

of N. with <|> , while the number K oc is the inner product of <*> with 

itself. Note that K g j is 1 x N where N is the number of Nj whose 

support intersects the support of (J> • For the above error estimates to be 

s 

uniformly valid as h+0, it is necessary that the support of <f> be fixed 

s 

independent of h. Thus N grows like 0(h - ^) as b*0. 

One of the early fears about the use of singular functions was that one 
might be giving up stability in order to gain greater accuracy. While it is 
true than the coefficient matrix in (4.9) has a far larger condition number 
that the standard stiffness matrix Kjj, this has not created ..problems in 

practice [25]. For example, in direct elimination the problems are isolated 
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in a one dimensional subspace. Indeed, the matrix in (4.9) admits the 
factorization 


(4.10) 


'll 




K 

ss 


L, , ^ 


V, , V 

11 


11 Is 

L . L 


V 

si ss. 


ss _ 


where 


(4.11) 


K 11 L 11 V 11 


(4.12) 


'Is 


L , . V. , K - L . V. , 
11 Is’ si si 11 


(4.13) 


K - L V, = L V . 
ss si Is ss ss 


Thus the bulk of the computation is in the factorization of the standard 

stiffness matrix Kjj in (4.11). The bordered parts L s j, V s j are obtained 

through the (stable) backsolves (4.12), and the only part of the calculation 

where the instability occurs in the lxl problem (4.13) (here L„„ = 1 

s s 

and V SS> K ss ■ L sl V ls are numbers). 

Also, note that if the stress intensity factor is the only 

quantity of interest, and this is frequently the case, then half of the 

backsolves can be omitted. Indeed, to get a, from (4.9) one backsolves 

n 


L lll 


= f 


L 03 
ss 


- L 


i* 


(4.14) 
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for y. and ui . Then 

(4.15) V a = u. 

ss 

The remaining backsolves, i.e., 

(4.16) V n + - y - V ls a, 
can be omitted. 

The rates of convergence for <t>^ and are the same as that 

observed for the case of grid refinement [12]. Thus asymptotically as h+0 
there is not a great deal of difference between the two approaches. Although, 
as noted earlier, the superiority of grid refinement emerges as adaptivity is 
incorporated into the approximation. 

In practical terms, the singular function approach has the advantage of 

fewer backsolves (in the case where only the intensity is desired), but 

suffers the disadvantage of requiring that not only inner products of nodal 

functions Nj be evaluated, but also inner products of the more complicated 

singular functions $ • Both of these effects contribute to the lower order 

terms in the overall work effect — considered as a function of — and thus 

ri 

both effects are negligible in the limit as h+0. 

On the other hand, for modest accuracy requirements (h not small) the 
savings in backsolves tends to be a significant effect. It is exactly in 
these cases where singular element methods have proven to be useful. To cite 
a specific numerical example, consider the problem described in Figure 4.1. 
This problem has been used as a test problem for a large number of different 
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nuraerical methods ([25] - [30]), and it is known that the stress intensity 

factor is 

o = .1917 

(to the number of decimals shown). Using a single singular function on a 

uniform grid with h « ^ (65 unknowns) gives 

o h = .1862 

while a 0-grid refinement with 16 unknowns (and over three times the CPU 

cycles and storage requirements) gave 

a, = .1621 
n 

with finer grids; however, the differences in overall work for a given 

accuracy rapidly disappears [27]. 

Another special case where singular elements have proven useful is what 
could be called an "h-p version" of the singular element method. Here one 
uses more singular functions than would be actually needed, not only to 

substract out the singularity but also to approximate. For example, in the 
problem cited above, the use of 8 singular functions and a uniform grid with 
h = -j- gives 

o, = .1916 

h 


(see [27] for this and other examples) 
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